{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'rasterio'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-10f54538d97a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mosgeo\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mgdal\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcrs\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mCRS\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwarp\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mcalculate_default_transform\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;32mfrom\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwindows\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mget_data_window\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'rasterio'"
     ]
    }
   ],
   "source": [
    "from osgeo import gdal\n",
    "import rasterio\n",
    "from rasterio.crs import CRS\n",
    "from rasterio.warp import calculate_default_transform\n",
    "from rasterio.windows import get_data_window\n",
    "from rasterio.vrt import WarpedVRT\n",
    "from rasterio.enums import Resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import cv2 as cv\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import corner_peaks, corner_harris"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RGBN_dir = \"E:\\\\datasets\\\\spacenet\\\\train\\\\MS\\\\SN6_Train_AOI_11_Rotterdam_MS_20190804111224_20190804111453_tile_8679.tif\"\n",
    "SAR_dir = \"E:\\\\datasets\\\\spacenet\\\\train\\\\SAR-Intensity\\\\SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190804111224_20190804111453_tile_8679.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RGBN_raster = rasterio.open(RGBN_dir)\n",
    "SAR_raster = rasterio.open(SAR_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RGBN_raster.crs.to_epsg()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_virtual_warped_raster(raster, projection=\"EPSG:3035\", resolution=2.5):\n",
    "        epsg_code = int(projection.split(':')[1])\n",
    "\n",
    "        # Broken PRISM images\n",
    "        if raster.crs is None:\n",
    "            raster._crs = CRS.from_epsg(epsg_code)\n",
    "\n",
    "        if raster.crs is not None and raster.crs.to_epsg() == epsg_code:\n",
    "            return raster\n",
    "        \n",
    "        dst_transform, dst_width, dst_height = calculate_default_transform(raster.crs, projection, raster.width, raster.height, *raster.bounds, resolution=(resolution, resolution))\n",
    "        vrt_opts = {\n",
    "            'resampling': Resampling.nearest,\n",
    "            'crs': projection,\n",
    "            'transform': dst_transform,\n",
    "            'height': dst_height,\n",
    "            'width': dst_width,\n",
    "            }\n",
    "\n",
    "        return WarpedVRT(raster, **vrt_opts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_from_geotif(src, band, roi=None):\n",
    "    '''\n",
    "    读取包括src在内的方形窗口（可能包含外侧的nodata值）\n",
    "    \n",
    "    Return:\n",
    "        img, (width, height), transform\n",
    "    '''\n",
    "    if roi is None:\n",
    "        roi = get_data_window(src.read(band, masked=True))\n",
    "\n",
    "    img = src.read(band, window=roi)\n",
    "    width = src.width\n",
    "    height = src.height\n",
    "    transform = src.transform\n",
    "\n",
    "    return img, (width, height), transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FakeCVFpt:\n",
    "    def __init__(self, xy):\n",
    "        self.size = 1\n",
    "        self.response = 1\n",
    "        self.pt = (xy[1], xy[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_keypoints(img, scheme=\"SURF\", radius=None):\n",
    "    \n",
    "    cnrs = corner_peaks(corner_harris(img), min_distance=radius)\n",
    "    kps = [FakeCVFpt(xy) for xy in cnrs]\n",
    "\n",
    "    return kps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RGBN_tif = create_virtual_warped_raster(RGBN_raster)\n",
    "RGBN_warped, _, _, _ = load_from_geotif(RGBN_tif)\n",
    "RGBN_raster = RGBN_raster.read()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "40d1469076e660c9bde21875d026b67d0c7c90423db199e68f2f3721b842c843"
  },
  "kernelspec": {
   "display_name": "Python (SOMatch)",
   "language": "python",
   "name": "somatch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
